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Abstract 


We report on work carried out to develop a 3-D coupled Finite Volume/BEM-based temperature 
forward/flux back (TFFB) coupling algorithm to solve the conjugate heat transfer (CHT) which 
arises naturally in analysis of systems exposed to a convective environment. Here, heat 
conduction within a structure is coupled to heat transfer to the external fluid which is convecting 
heat into or out of the solid structure. There are two basic approaches to solving coupled fluid- 
structural systems. The first is a direct coupling where the solution of the different fields is 
solved simultaneously in one large set of equations. The second approach is a loose coupling 
strategy where each set of field equations is solved to provide boundary conditions for the other. 
The equations are solved in turn until an iterated convergence criterion is met at the fluid-solid 
interface. The loose coupling strategy is particularly attractive when coupling auxiliary field 
equations to computational fluid dynamics codes. We adopt the latter method in which the BEM 
is used to solve heat conduction inside a structure which is exposed to a convective field which 
in turn is resolved by solving the NASA Glenn compressible Navier-Stokes finite volume code 
Glenn-HT. The BEM code features constant and bi-linear discontinuous elements and an ILU- 
preconditioned GMRES iterative solver for the resulting non-symmetric algebraic set arising in 
the conduction solution. Interface of flux and temperature is enforced at the solid/fluid interface, 
and a radial-basis function scheme is used to interpolated information between the CFD and 
BEM surface girds. Additionally, relaxation is implemented in passing the fluxes from the 
condution solution to the fluid solution. Results from a simple test example are reported. 


1 Introduction 


As field solvers have matured, coupled field analysis has received much attention in an effort to 
obtain computational models ever-more faithful to the physics being modeled, see for instance 
[1]. The coupled field problem which we address is the conjugate heat transfer (CHT) problem 
arising commonly in practice: time dependent or time independent convective heat transfer over 
coupled to conduction heat transfer within a solid body, see Fig. (1). Examples of CHT include 
analysis of automotive engine blocks, fuel ejectors, cooled turbine blade/vanes, nozzle or 
combustor walls, or thermal protection system for re-entry vehicles. Applications of interest are 
then any thermal system in which multi-mode convective/conduction heat transfer is of particular 
importance to thermal design, and thus CHT arises naturally in most instances where external 
and internal temperature fields are coupled. 



Figure 1 : CHT problem: external convective heat transfer coupled to heat conduction 

within the solid. 

However, conjugacy is often ignored in numerical simulations. For instance, in analysis of 
turbomachinery, separate flow and thermal analyses are typically performed, and, a constant wall 
temperature or heat flux boundary condition is typically imposed for the flow solver. Convective 
heat transfer coefficients are then obtained from the flow solution, and these are provided to the 
conduction solver to determine the temperature field and eventually to perform further thermal 
stress analysis. The shortcomings of this approach which neglects the effects of the wall 
temperature distribution on the development of the thermal boundary layer are overcome by a 
CHT analysis in which the coupled nature of the field problem is explicitly taken into account in 
the analysis. 

There are two basic approaches to solving coupled field problems. In the first approach, a direct 
coupling is implemented in which different fields are solved simultaneously in one large set of 
equations. Direct coupling is mostly applicable for problems where time accuracy is critical, for 
instance, in aero-elasticity applications where the time scale of the fluid motion is on the same 
order as the structural modal frequency. However, this approach suffers from a major 
disadvantage due to the mismatch in the structure of the coefficient matrices arising from 
BEM, FEM and/or FVM solvers. That is, given the fully populated nature of the BEM 


coefficient matrix, the direct coupling approach would severely degrade numerical efficiency of 
the solution by directly incorporating the fully populated BEM equations into the sparsely banded 
FEM or FVM equations. A second approach which may be followed is a loose coupling strategy 
where each set of field equations is solved separately to produce boundary conditions for the 
other. The equations are solved in turn until an iterated convergence criterion, namely continuity 
of temperature and heat flux, is met at the fluid-solid interface. The loose coupling strategy is 
particularly attractive when coupling auxiliary field equations to computational fluid dynamics 
codes as the structure of neither solver interferes in the solution process. 

There are several algorithmic approaches which can be taken to solve the coupled field problems. 
Most commonly used methods are based on either finite elements (FEM), finite volume methods 
(FVM), or a combination of these two field solvers. Examples of such loosely coupled 
approaches applied to a variety of CHT problems ranging from engine block models to 
turbomachinery can be found in Comini et al. [2], Shyy and Burke [3], Patankar [4], Kao and 
Liou [5], Hahn et al. [6], Bohn et al. [7], and in Tayla et al. [8] where multi-disciplinary 
optimization is considered for CHT modeled turbine airfoil designs. Hassan et al. [9] develop a 
conjugate algorithm which loosely couples a FVM-based Hypersonic CFD code to an FEM heat 
conduction solver in an effort to predict ablation profiles in hypersonic re-entry vehicles. Here, 
the structured grid of the flow solver is interfaced with the un-structured grid of heat conduction 
solvers in a quasi-transient CHT solution tracing the re-entry vehicle trajectory. Issues in loosely 
coupled analysis of the elastic response of solid structures perturbed by external flowfields 
arising in aero-elastic problems can be found in Brown [10]. In either case, the coupled field 
solution requires complete meshing of both fluid and solid regions while enforcing solid/fluid 
interface continuity of fluxes and temperatures, in the case of CHT analysis, or displacement and 
traction, in the case of aero-elasticity analysis. 

A different approach taken by Li and Kassab [11,12] and Ye et al. [12] who develop a HCT 
algorithm which avoids meshing of the solid region to resolve the heat conduction problem. In 
particular, the method couples the boundary element method (BEM) to an FEM Navier-Stokes 
solver to solve a steady state compressible subsonic CHT problem over cooled and uncooled 
turbine blades. Due to the boundary-only discretization nature of the BEM, the onerous task of 
grid generation within intricate solid regions is avoided. Here, the boundary discretization 
utilized to generate the computational grid for the external flow-field provides the boundary 
discretization required for the boundary element method. In cases where the solid is multiply- 
connected, such as a cooled turbine blade, the interior boundary surfaces must also be 
discretized; however, this poses little additional effort. Moreover, in addition to eliminating 
meshing the solid region, this BEM/FVM method offers an additional advantage in solving CHT 
problems which arises from the fact that nodal unknowns which appear in the BEM are the 


surface temperatures and heat fluxes. Consequently, solid/fluid interfacial heat fluxes which are 
required to enforce continuity in CHT problems are naturally provided by the BEM conduction 
analysis. This is in contrast to the domain meshing methods such as FVM and FEM where heat 
fluxes are computed in a post-processing stage by numerical differentiation. Efe et al. [14,15] 
adopted this approach in further studies of CHT in incompressible flow in ducts subjected to 
constant wall temperature and constant heat flux boundary conditions. Kontinos [16] also 
adopted the BEM/FVM coupling algorithm to solve the CHT over metallic thermal protection 
panels at the leading edge of the X-33 in a Mach 15 hypersonic flow regime. Rahaim and Kassab 
[17] and Rahaim et al. [18] adopt a BEM/FVM strategy to solve time-accurate CHT problems for 
supersonic compressible flow, and they present experimental validation of this CHT solver. In 
their studies, the dual reciprocity BEM [19] is used to for transient heat conduction, while a cell- 
centered FVM is chosen to resolve the compressible turbulent Navier-Stokes equations. 

We now present work carried out under this grant which extends the loosely coupled BEM/FVM 
approach to solving CHT problems in steady state for 3-D problems. Here, we couple a 3-D 
BEM conduction solver with the NASA Glenn multi-block FVM Navier-Stokes convective heat 
transfer code, Glenn-HT. This density based FVM code is a robust and well-proven code 
developed specifically for fluid-flow and heat transfer analysis in turbomachinery applications. 
The methodology adopted in this work is given as well as the information passing process. 
Results are presented for a test-case configuration to be used in future laboratory experiments 
which will serve as experimental validation of the CHT solver. 

2 Governing Equations 

The governing equations for the mixed field problem under consideration are reviewed. The 
CHT problems arising in turbomachinery involves external flow fields that are generally 
compressible and turbulent, and these are governed by the compressible Navier-Stokes equations 
supplemented by a turbulence model. Heat transfer within the blade is governed by the heat 
conduction equation, and linear as well as non-linear options are considered. However, fluid 
flows within internal structures to the blade, such as film cooling holes and channels, are usually 
low-speed and incompressible. Consequently, density-based compressible codes are not directly 
applicable to modeling such flows, unless low Mach number pre-conditioning is implemented, 
see Turkel [25,26]. It is also noted that the Glenn-HT code is specialized to turbomachinery 
applications for which air is the working fluid and which is modeled as an ideal gas. 

2.1 The Flow Field 

The governing equations for the external fluid flow and heat transfer are the compressible 
Navier-Stokes equations, which describe the conservation of mass, momentum and energy. These 
can be written in integral form as 


n 


r 


n 


(i) 


where Q denotes the volume, V denotes the surface bounded by the volume ft, and n is the 
outward-drawn normal. The conserved variables are contained in the vector W = (p, pu, pv y 
pw , pe, pk , pcu) , where, p, u, u;, e, k y u) are the density, the velocity components in x- y y -, 
and ^-directions, and the specific total energy, k and uj are the kinetic energy of turbulent 
fluctuations and the specific dissipation rate of the two equation k-u Wilcox turbulence model 
[23, 24] is adopted in Glenn-HT. The vectors £ and J are convective and diffusive fluxes 
respectively, g is a vector containing all terms arising from the use of a non-inertial reference 
frame as well as production and dissipation of turbulent quantities. The fluid is modeled as an 
ideal gas. A rotating frame of reference can be adopted for rotating flow studies. The effective 
viscosity is given by 

P = P/ + Pi (2) 

where p t ~ pk/u. The thermal conductivity of the fluid is then computed by a Prandtl number 
analogy where 


7 P* Pi 

7 — 1 Pri Pr t 


( 3 ) 


and Pr is the Prandtl number and 7 is the specific heat ratio. The subscribes l and t refers to 
laminar and turbulent values respectively. 


2.2 The Heat Conduction Field 

In the CHT model, the NS equations are solved to steady state by a time marching scheme. As 
physically realistic time-dependent solutions are not sought, a quasi-steady heat conduction 
analysis using the BEM is performed at a given time level. As such, only steady-state heat 
conduction is considered. The governing equation is 

V ■ [k(T s )VT s ] = 0 (4) 

where, T 6 denotes the temperature of the solid, and k s is the thermal conductivity of the solid 
material. If the thermal conductivity is taken as constant, then the above reduces to the following 
equation for the temperature 

fcV 2 T = 0 (5) 

When the thermal conductivity variation with temperature is an important concern, the 
nonlinearity of the heat conduction equation can readily be removed by introducing the classical 
Kirchhoff transform, f7(T), which is defined as 


( 6 ) 
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where T 0 is the reference temperature and k Q is the reference thermal conductivity. The 
transform and its inverse are readily evaluated, either analytically or numerically, and the heat 
conduction equation transforms to a Laplace equation for the transform parameter U(T), that is 
the nonlinear heat conduction problem transforms to 

k 0 V 2 U = 0 (7) 

and the boundary conditions of the first and second kind transform linearly as 

U — U{T S ) and k a ^- = k s ^ (8) 


Thus, the nonlinear problem in T subjected to any combination of 1st or 2nd kind boundary 
conditions transforms to a linear problem in the transform parameter U. The heat conduction 
equation thus reduces to the Laplace equation in any case, and this equation is readily solved by 
the BEM. 

In the conjugate problem, continuity of temperature and heat flux at the blade surface, T, 
must be satisfied: 


Tf = Is 

kf d J_ i=- k ?l± 

Lf dn ks On 


( 9 ) 


Here, Tj is the temperature computed from the N-S solution, T s is the temperature within the 
solid which is computed from the BEM solution, and d/dn denotes the normal derivative. Both 
first kind and second kind boundary conditions transform linearly in the case of temperature 
dependent conductivity. As will be explained later, in such a case, the fluid temperature is used to 
evaluate the Kirchhoff transform and this is used a boundary condition of the first kind for the 
BEM conduction solution in the solid. Subsequently the computed heat flux, in terms of U, is 
scaled to provide the heat flux which is in turn used as an input boundary condition for the flow- 
field. 


3 Field Solver Solution Algorithms 

A brief description of the Glenn-HT code is given in this section. Details of the code and its 
verification can be found in [21,22]. The heat conduction equation is solved using BEM, and 
details are provided for this conduction solver. 

3.1 Navier-Stokes Solver 

A cell-centered FVM is used to discretize the NS equations. Eq. (1) is integrated over a 
hexahedral computational cell with the nodal unknowns located at the cell center k). The 



convective flux vector is discretized by a central difference supplemented by artificial dissipation 
as described in Jameson et al. [27]. The artificial dissipation is a blend of first and third order 
differences with the third order term active everywhere except at shocks and locations of strong 
pressure gradients. The viscous terms are evaluated using central differences. The resulting finite 
volume equations can be written as 
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( 10 ) 


where W ... is the cell-volume averaged vector of conserved variables, q and d ... are the 
net flux and dissipation for the finite volume obtained by surface integration of Eq. (1), and^s 
is the net finite source term. The above is solved using a time marching scheme based on a fourth 
order explicit Runge-Kutta time stepping algorithm. The steady-state solution is sought by 
marching in time until the dependent variables reach their steady-state values, and, as such, 
intermediate temporal solutions are not physically meaningful. In this mode of solving the 
steady-state problem, time-marching can be viewed as a relaxation scheme, and local time- 
stepping and implicit residual smoothing are used to accelerate convergence. A multi-grid option 
is available in the code, although it was not used in the results reported herein. 


3.2 Heat Conduction Boundary Element Solution 

The heat conduction equation reduces to the same governing Laplace equation in the temperature 
or the Kirchhoff transform. In the boundary element method this governing partial differential 
equation is converted into a boundary integral equation (BIE), see[28-31], as 

C(0T(0+ I T(x)q*(x,OdS(x)= j q(x)T*(x,OdS(x ) (1 1) 

Js Js 

where 5(x) is the surface bounding the domain of interest, £ is the source point, x is the field 
point, q{x) = — kdT/dn is the heat flux, T*(x, £) is the so-called fundamental solution, and 
g*(x,£) is its normal derivative with d/dn denoting the normal derivative with respect to the 
outward-drawn normal. The fundamental solution (or Green free space solution) is the response 
of the adjoint governing differential operator at any field point xdue to a perturbation of a Dirac 
delta function acting at the source point <£. In our case, since the steady state heat conduction 
equation is self-adjoint, we have 

fcV 2 T*(x,£)= (12) 

Solution to this equation can be found by several means, see for instance Liggett and Liu [32], 
Morse and Feschbach [33], and Kellog [34], as 


(13) 
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where r(x, <Q is the Euclidean distance from the source point £. The free term C(£) can be shown 
analytically to be 

C(0=<f -k dT ^ X '° dS(x) ( 14 ) 

JS(x ) an 

Moreover, introducing the definition of the fundamental solution in the above, it can be readily 
be determined that C(£) is the internal angle (in degrees in 2-D and in steradians in 3-D) 
subtended at source point divided by 27rin 2-D and by 47r in 3-D when the source point £is on 
the boundary and takes on a value of one when the source point £ is at the interior. Consequently, 
the free term takes on values 1 > C (0 > o. 

In the BEM, the BIE is discretized using two levels of discretization; 

1 . the surface 5 is discretized into a series of e = 1,2 ...N elements. This is traditionally 
accomplished using polynomial interpolation, bilinear and quadratic being the most 
common. In general, 

N 

s = J2 Ase ( 15 ) 

e— 1 

and on each surface element S e the geometry is discretized using local shape functions 
Q in terms a homogenous coordinates (77, () which each take on values between 
[-1,1] as 

NGE 

0 = 2 N * fa O x fc ( 16 ) 
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Here, (x^, y*, z^) denote the location of the k = 1, 2 ...NGE boundary nodes used to 
define the element geometry. 

2. the distribution of the temperature and heat flux is modeled on the surface. This is 
usually accomplished using polynomial interpoaltion as well. Common discretizations 


include: constant (where the mean value of T and q) are taken on an element surface, 
blinear, or bi-quadratic. In general, 

NPE 

T'( V , C) = £>/(»?, 0*/ (17) 

J=1 

NPE 

q e (v,0 = ^2 N j{vX)Q^ 

j= 1 

It is noted that the order of discretization of the temperature and heat flux need not be 
the same as that used for the geometry, leading to sub-parametric (lower order than that 
used for the geometry), iso-parametric (lower order than that used for the geometry), 
and super-parametric (lower order than that used for the geometry) discretizations. 
Moreover, the temperature and heat flux are discretized using j — 1,2, ...NPE 
discrete nodal values whose location within the element e can be chosen to 


(a) coincide with the location of the geometric nodes: continuous elements. 

(b) be located offset from the geometric nodes: discontinuous elements. 

We choose to employ bi-linear discontinuous elements as they provide high levels of accuracy in 
computed heat flux values especially at sharp comers regions where first kind boundary 
conditions are imposed without resorting to special treatment of comer points required by 
continuous elements[30]. Such comer regions are often encountered in industrial problems and 
first kind boundary conditions are imposed there in the CHT algorithm to be shortly described. 
Moreover, the use of discontiuous elements throughout the BEM model eliminates much of tyhe 
overhead associated with continuous elements, in particular, there is no need to generate, store, 
or access a connectivity matrix when using discontinuous elements. Details of the discretization 
employed in the BEM code are provided in the Appendix. 


Following standard BEM discretization the BIE leads to the following 


N NPE 
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The surface integrals in the above equation depend purely on the local geometry of the element 
and the location ofthe source point £. These are evaluated numerically using Gauss quadratures. 
Upon collocation of the above at every boundary node where the temperature and heat flux are 
define, the following algebraic form is obtained: 


[tf]{T s } = [G}{ q ,} 


(19) 


Here the influence matrices [H] and [G] are evaluated numerically using quadratures. Details of 
the numerics are provided in the Appendix. These are solved subject to either of the following 
boundary conditions: 

(a) at the external bounding wall: T s |p e — T f (20) 

(b) at internal cooling hole surfaces: — fc s |^| re — h[T s — ]| Fc (21) 

Here, T is the wall temperature computed from the N-S solution, T s is the wall temperature 
computed from the BEM conduction solution, the external boundary is denoted by r c , while F c 
denotes the convective internal cooling hole boundaries. At this point the Glenn-HT code does 
not have low Mach number capabilities. As such, modeling of convection in the cooling passages 
must be accomplished using correlations and not simulated using the NS solver (although this is 
perfectly possible and would be consistent with the CHT philosophy with proper extensions of 
the NS solver). In addition, in certain cases, a full CHT model is not practical in the design stage. 
For instance, modeling of cooling passages in a turbine blade with intricate cooling schemes such 
as jet impingement cooling or turbulence enhancing ribs pose a serious computational challenge, 
and these are often better modeled using correlation equations. Consequently, a partial CHT 
model can be carried out with internal bounding walls modelled using the standard convective 
boundary condition at internal cooling hole surfaces. 

The advantage of the BEM formulation over finite difference or finite element formulations is 
that no interior mesh is generated and the surface heat flux is computed in the solution. A 
GMRES iterative solver with an ILU pre-conditioning is used to solve the BEM equations. 

3.3 CHT Algorithm 

The Navier-Stokes equations for the external fluid flow and the heat conduction equation for heat 
conduction within the solid are interactively solved to steady state through a time-marching 
algorithm. The surface temperature obtained from the solution of the Navier-Stokes equations is 
used as the boundary condition of the boundary element method for the calculation of heat flux 
through the solid surface. This heat flux is in turn used as a boundary condition for the Navier- 
Stokes equations in the next time step. This procedure is repeated until a steady-state solution is 
obtained. In practice, the BEM is solved every few cycles of the FVM, say every two hundred, to 
update the boundary conditions, as intermediate solutions are not physical in this scheme. This is 
referred to as the temperature forward/flux back (TFFB) coupling algorithm as outlined below: 


• FVM Navier-Stokes solver: 

1 . Begins with initial adiabatic boundary condition at solid surface. 

2. Solves compressible NS for fluid region. 

3. Provides temperature distribution to BEM conduction solver 
after a number of iterations. 

4. Receives flux boundary condition from the BEM as input for 
next set of iterations. 


• BEM conduction solver: 

1. Receives temperature distribution from FVM solver. 

2. Solves steady-state conduction problem. 

3. Provides flux distribution to FVM solver. 


The transfer of heat flux from the BEM to the FVM solver is accomplished after under- 
relaxation. 
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= P Qold M + (1 - 
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( 22 ) 


with P taken as 0.2 in all reported calculations. The choice of the relaxation parameter is through 
trial and error. In certain cases, it has been our experience that a choice of larger relaxation 
parameter can lead to non-convergent solutions[35]. The process is continued until the NS solver 
converges and wall temperatures and heat fluxes converge, that is until Eq. (23) is satisfied 
within a set tolerance 

111,-1,11 <^r (23) 


where and e q are taken as 0.001. 


It should be noted that alternatively, the flux may be specified as a boundary condition for the 
BEM code leading to a flux forward temperature back (FFTB) approach. However, when a 
fully conjugate solution is undertaken, this would amount to specifying second kind boundary 
conditions completely around the surface of a domain governed by an elliptic equation, 
resulting in a non-unique solution. Thus, the TFFB algorithm avoids such a situation. A blend 
of the TFFB and FFTB may resolve this problem and is the subject of ongoing work. Results 
of these studies will be reported elsewhere. 


3.4 Interpolation Between BEM and FVM Grids 

A major issue in information transfer between CFD and BEM is the difference in the levels of 
discretization between the two meshes employed in a typical CHT simulation. Accurate 


resolution of the boundary layer requires a FVM surface grid which is much too fine to be 
used directly in the BEM. A much coarser surface grid is typically generated for the BEM 
solution of the conduction problem. The disparity in meshes is illustrated for the test section 
considered in this paper, see Fig. (2): the FVM grid in Fig. (3) used for the external flow 
solution is obviously much finer than the coarser BEM grid in Fig. (4) which was used for 
internal conduction analysis. The disparity between the two grids requires a general 
interpolation of the surface temperature and heat flux between the two solvers as it is not 
possible in general to isolate a single BEM nodes and identify a set of nearest FVM nodes. 
Indeed in certain regions where the CFD mesh is very fine, a BEM node can readily be 
surrounded by tens or more FVM nodes. 

Radial basis function (RBF) interpolation [19,20] can naturally be adopted for this purpose. 
Consider Fig. (5), here the location of a BEM node is identified on the right-hand-side by a star- 
like symbol. Let the position of the BEM node of interest be denoted by r, then, the value of the 
temperature at that location is interpolated using RBF’s with poles located at the position r* of 
each of the FVM surface nodes lying within a sphere of radius Rm ax centered about r 

NPS 

T (?) =^2 a ifi( r ) ( 24 ) 

1=1 

where, NPS is the total number of FVM nodes contained within the sphere, r — |r — r*| 
denotes the Euclidean distance, and /*(r) are radial basis functions. Here, we use the standard 
conic function 

f(r) = l+r (25) 

Collocating Eq. (24) at all i = 1,2. ..NPS nodes of the finite volume mesh surrounding the 
BEM node of interest, the expansion coefficients are found by solving the equation 

T = Fa (26) 

where T is the vector of FVM computed wall temperatures, and the interpolant matrix jf is 
known from RBF theory to be well conditioned and readily inverted. In all calculations, the 
maximum radius R m ax of the sphere is set to 5% of the maximum distance within the solid 
region. This limit may be adjusted to suit the problem at hand. 

4 Numerical Results 

We now report on a preliminary simulation used to verify our algorithm. A 3-D model is made of 
a 2-D configuration used in an experiment set up to simulate heat transfer conditions in a cooled 
turbine blade tip and investigate the importance of conjugacy, see Fig. 2. Heated air at 319. 5K 
enters horizontally at the left end of the channel, flows over the simulated tip gap, and out the 
bottom of the channel. The block is made of stainless-steel (k ~ 14.9 W/mK, 


p — 8.03x10 3 fcp/m 3 , c = 502.48 J/kgK) and is cooled by a laminar flow of water at 286 K. 
The film coefficient in the cooling channels is calculated as 536 W/m 2 K. In this simulation, 
convective boundary conditions are used to model flow in the channel as a full external flow and 
internal flow CHT solution was not carried out, due to the fact that the CFD code is specialized 
for turbomachinery applications and can only model air as a working fluid. All walls not 
exposed to the flow are adiabatic. The height of inlet channel is 0.02m, the height of the outlet 
channel is 0.0025m, width of passage between block and wall is 0.005m, the length of inlet 
channel up to block is 0.40315m, the length of outlet channel measured from block to exit is 
0.40615m. Total pressure at inlet is atmospheric and the back pressure at exit is 0.92 x inlet total 
pressure (i.e., p/P_inlet = 0.92). A 3-D model was constructed to model the centerline of the 
block. Four finite volume cells were used to model the width and the surface grid at the block are 
shown in Figs. 3 and 4. The FVM mesh uses 1 104 surface cells, with a total of 200,000 finite 
volume cells. The BEM surface grid used 946 uniformly spaced boundary elements. The CHT 
solution was run to reduce residuals in density to 1 .0xl0~ 6 and in energy to 1 .OxlO -9 . Results for 
the CHT predicted block surface temperature, flow passage temperature and Mach number are 
displayed in Figs. 6 and 7. The Mach number in center of inlet channel is approximately 0.06, 
and typical Mach number in outlet channel is approximately 0.068 (variation over the height of 
the channel due to the flow over the block). 

5 Conclusions 

A combined FVM/BEM method has been developed to solve the conjugate problem in CHT 
analysis. As a boundary only grid is used by the BEM, the computational time for the heat 
conduction analysis is insignificant compared to the time used for the NS analysis. The proposed 
method produces realistic results without using arbitrary assumptions for the thermal condition at 
the conductor surface. In practice, turbomachinery components such as modem cooled turbine 
blades which often contain upwards of five hundred film cooling holes and intricate internal 
serpentine cooling passages with complex convective enhancement configurations such as 
turbulating trip strips, pose a real computational challenge to BEM modeling. It is proposed to 
extend the current work by implementing either fast multipole-accelerated BEM or adoption of 
iterative block solvers for the BEM in order to address large-scale problem. 
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Figure 2: Cross-section of experiment set-up which is also used in numerical simulation. 

Upstream channel extends 10 hydraulic diameters upstream of the block. 



Figure 3; Finite Volume surface mesh: cell centered finite volumes with four cells in 
z-direction and total of 1 104 surface cells. 


Figure 4: 


Surface BEM discretization for the block with 946 equally spaced bilinear 
elements distributed over the surface of the block. 




Figure 5: Transfer of nodal values from FVM and BEM (and back) independent surface 

meshes is performed with a compactly supported radial-basis-function 
interpolation. 



Figure 6: Surface temperature distribution from the converged conjugate solution 

(temperature distribution plotted from the BEM solution). 


Figure 7: 
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(C) 


CHT predicted (a) flow passage temperatures and Mach numbers, (b) top comer 


Mach number, and (c) top comer temperatures. 


APPENDIX 


3D BEM FORMULATION FOR STEADY-STATE HEAT CONDUCTION 

The point of departure is the BIE for the steady state heat conduction equation, 

0(0^)+ I T(x)q* (x, £)dS(x)= I q(x)T*(x,OdS(x ) (A.l) 

Js Js 

where x denotes the 3-D coordinates, S(x) is the surface bounding the domain of interest, £ is 
the source point, x is the field point, q{x) = — kdT/dn is the heat flux, T*(x,^) is the so- 
called fundamental solution, and q*{x , £) is its normal derivative with d/dn denoting the normal 
derivative with respect to the outward-drawn normal. The fundamental solution (or Green free 
space solution) is the response of the adjoint governing differential operator at any field point x 
due to a perturbation of a Dirac delta function acting at the source point £. In our case, since the 
steady state heat conduction equation is self-adjoint, we have 

w 2 r(x,o= -s(x,o (a.2) 

Solution to this equation can be found by several means which yield 

T *OuO = - 7TT lnr (' GO in2 ‘ D (A.3) 


47 rk r(x , £) 

where r(x, £) = \x — £| is the Euclidean distance from the source point £ to the field point x. The 
free term C(£) can be shown analytically to be 

c(d)=<£ -k dT '^ dS (A.4) 

Js{x ) <vn 

Moreover, introducing the definition of the fundamental solution in the above, it can be readily 
be determined that C(£) is the internal angle (in degrees in 2-D and in steradians in 3-D) 
subtended at source point divided by 27rin 2-D and by 47T in 3-D when the source point £ is on 
the boundary and takes on a value of one when the source point £ is at the interior. 
Consequently, the free term takes on values 1 > C(£) > 0. 


In the BEM, the BIE is discretized using two levels of discretization: 


1. the surface S is discretized into a series of e = 1,2 ...N elements. This is 

traditionallyaccomplished using polynomial interpolation, bilinear and quadratic 
being the most common. In general, 

S = ^AS e (A. 5) 

e—\ 

and on each surface element A S e {x) the geometry is discretized using local shape 
functions A 7 ^ ( 77 , £) in terms a homogenous coordinates ( 77 , £) which each take on values 
between [ — 1 , 1 ] as 

NGE 

x e {' n ,o = Y. N ^Q x k (A. 6 ) 

1 

NGE 

y e M = Y^ N k( r hOyk 

1 

NGE 

2 ! 0j,o = X> t e 07.c)* 

k = 1 

Here, (xk,yk, Zk) denote the location of the k = 1,2 ...NGE boundary nodes used 
to define the element geometry. 


2 . the distribution of the temperature and heat flux is modeled on the surface. This is 
usually accomplished using polynomial interpoaltion as well. Common 
discretizations include: constant (where the mean value of T and q) are taken on an 
element surface, blinear, or bi-quadratic. In general, 

NE 

T e (r 7,0 = ^1^(77,077 (A.7) 

j=l 

NE 

q e {y,Q = 

j = 1 

It is noted that the order of discretization of the temperature and heat flux need not 
be the same as that used for the geometry, leading to sub-parametric (lower order 
than that used for the geometry) , iso-parametric (lower order than that used 

for the geometry), and super-parametric (lower order than that used for the 
geometry) discretizations. Moreover, the temperature and heat flux are 


discretized using j = 1,2, ...NE discrete nodal values whose location within the 
element e can be chosen to 

(a) coincide with the location of the geometric nodes: continuous elements. 

(b) be located offset from the geometric nodes: discontinuous elements. 

Tw 7 o types of discontinuous boundary elements are available in the BEM code used in the 
project: constant elements (sub-parameteric) and bilinear isoparametric. All CHT calculations are 
reported using bi-linear discontinuous elements. 

Constant Element 


In the constant element, which is a sub-parametric element, the field variables, T and q, are 
modeled as constant across each element while the geometry is represented locally as bi-linear 
planes. Figure A.l below shows a typical constant boundary element along with its transformed 
representation in the local 77 — £ coordinate system. 



Figure A. 1 . Bi-linear subparametric boundary element. 


Notice that the geometric nodal locations of the element are ordered counterclockwise such that 
the normal vector always points outwards from the domain of the problem. The global coordinate 
system (z, y, z) is transformed into a local coordinate system (rj, Q using the following bi-linear 
shape functions relationships as, 


(A.8) 


z e (^o = c)xk 

k = 1 

y e {v, 0 = ^2^k(v, 0 Vk 

k=l 

z e {y,0 = 

i t=i 

where the four bi-Iinear shape functions are defined as follows, 

^(^0 = ^(1 -^7)(1 -0 (A.9) 

N 2 e (v,0 = ^(1 + T7)(1-C) 

N, e (v, C)« j(l + »/)(! + 0 
N 4 e (vX) = j(i - ??)(i + 0 

The temperature and heat flux are modeled as constant with the node located at the geometric 
center of the boundary element, thus 

T e (t7,C) = Tf and q e ( r } ,Q=q e j (A. 10) 

Clearly, the temperature and heat flux are discontinuous at the element interfaces. Thus, a 
constant element is termed discontinuous. 

Introducing the above discretization in the BIE in Eq. (A.l), noting that NPE = 1, and 
collocating the discretized BIE at each of the boundary nodes & there results 

N N 

cm rm + = ypm (A.10 

j= 1 3 = 1 

where the influence coefficients H tj and G XJ are defined as 

Hij= I I q*{x, £i) dS (A. 12) 

J Jas, 

Gij = I I T*(x, £,•) dS 
J JaSj 

These are evaluated numerically using Gauss-Legendre quadratures with an adaptive scheme to 
be discussed shortly. Although very simple in implementation, use of the the above constant 
element formulation does not lead to satisfactory results in many cases. 


Bi-Linear Isoparametric Discontinuous Elements 


In this type of boundary element which is used in all CHT calculations, the field variables T and 
q are modeled with discontinuous bi-linear shape functions across each element while the 
geometry is represented locally as continuous bi-linear surfaces. Figure A. 2 below shows a 
typical bi-linear subparametric boundary element along with its transformed representation in the 
local coordinate system. 



Figure A. 2. Bi-linear isoparameteric discontinuous boundary element. 


Again, the global coordinate system ( x , y, z) is transformed into a local coordinate system (77, Q 
where the four bi-linear shape functions defined in Eq. (A.9). The field variables, T and q, are 
modeled to vary bi-linearly across the boundary element through the use of four discontinuous 
shape functions with nodes located at an off-set position of 12.5% from the edges of the element. 
The field variables and shape functions are described as follows: 

n>j,() = E^Oj; e (a. 1 3) 

j= 1 

Q e (v, 0 = 0 Qj 

j = 1 


with the bilinear shape functions M® defined as, 


(A. 1 4) 


Introducing the above discretization in the BIE in Eq. (A.l), and collocating the discretized BIE 
at each of the boundary nodes & there results 

N NPE N NPE 

c«,) m) + H i T > = £ £ G h i‘ 

e=l j=l j=l j =1 

where the influence coefficients H \ L and G ^ are defined as 

H ij = I I q*(.x,£i)M;(ri,QdS 

J J A5j 

G ij = J j^T(x^i)Mf(ri,OdS 

These influence coefficients are again evaluated numerically via quadratures 

Numerical Evaluation of the Influence Coefficients 

The process is illustrated by considering constant elements. Introducing the definition of bilinear 
representation of the geometry Eq. (A. 8) into the constant element influence coefficient 
definition, Eq. (A. 12), the influence coefficient integrals are explicitly, 

e y = /’ (A.i7) 

H,J= f J T\x(n,0,^\JAnX)\dr,di 

The Jacobian of the transformation over the j-th element, | Jj(g, £)| is 

Jj(rh 0 = \] 9w9X ~ ( A - 18 ) 

where g VTI , , and g v ^ are the components of the metric tensor defined as. 


(A.l 5) 
(A. 16) 


(A. 19) 


The metrics ||, |*r, and || are readily found by differentiation of x e (t], £), 

y e {rj , £), and z e (r), £) in Eq. (A. 8). Introducing the metrics into the Jacobian and simplifying 
leads to, 


M 1 7.0 = 


l 


<92 
dr] < 9 £ 


dz d y . 

dy dQ J 


, dz dx 

~ ~ — I + [jhjd£ 


dx dz . 2 , dx dy 

drj d( ) + \ dr] d(, 


dy dx , 2 
dr] d() 


(A. 20) 


In the expression for q* (x, &) there arises the need to evaluate the outward-drawn normal n, as 
by definition 

dT*{x,& 


<7*0, 6) = 




= 0, &) 


The components of the unit vector on each elements can be easily computed using, 

dy dz dz dy 


n T — 


n„ 


n~ — 


i dz dx dx dz . . . . 

= Us?"a7s?J / 


(A. 21) 


(A. 22) 


The numerical integration process necessary to obtain the influence coefficients in Eq. (A. 17) is 
performed by double Gaussian quadratures (Gauss-Legendre specifically) simultaneously along 
the two local axis rj and £, leading to the following form, 

NG NG 

Ga = Yj Y w - WmE ^ nm 5 C nmi j{Vnmi C nm )l (A.23) 

7i=l m— 1 
NG NG 

Hij = ^ ^ ^ ^ JT/ri Fj ( 77um ^ Cnm > £ i) I‘^j( 7 7nm5 Crim)| 

72— 1 771=1 


where NG is the number of Gaussian points or order of integration employed, r\ nm and ( nm are 
the locations of the Gaussian quadrature points (zeroes of the appropriate Legendre polynomials), 
and W n and W m are the quadrature weights [1,2]. 

The number of Gaussian points employed can be addapted to every integral depending of the 
variability of the integrand. The influence coefficients G are inversly proportional to the 
Euclidean distance between the field or integration point x and the collocation point and the 
influence factors H are inversly proportional to the square of the Euclidean distance between the 
field or integration point and the collocation point Therefore, as the collocation point is 
positioned closer to the integration element, the variability of the integrand increases requiring an 
increase in the number of Gaussian points, NG, for the integral approximation to provide a 
similar level of accuracy. Hence, a simple distance rule can be heuristically employed to change 
the number of Gaussian points depending on how far the collocation point is to the integration 
element. 

However, simply increasing the number of Gaussian integration points is in some cases not 
enough to obtain an accurate approximation to the integral. This is precisly the case when the 
collocation point is extremely close to the integration element. In such cases a 
subsegmentation of the element is required and can be performed in a similar fashion as to the 
increase of Gaussian points, that is, the closer the collocation point is to the element, the more 
subdivisions are made to the element with a fixed number of Gaussian points for each 
subelement. 

The particular case in which the collocation point is located over the integration element must 
be treated with caution. For this case the integrand becomes singular lacking an accurate integral 
approximation through a regular Gaussian rule. Even though the integral for the influence factor 
Ha is strongly singular because its integrand is inversely proportional to the distance squared 
(1/r 2 ), this integral need not be computed directly but instead using the equipotential relation by 
suming the off-diagonal terms[3-5], that is 

N 

H»= H H (A.24) 

The integral for the influence coefficient G it is however weakly singular because its integrand is 
inversely proportional to the distance (1/r), therefore the singularity may be avoided through an 
additional transformation of the local coordinate system. Figure A. 3 shows the primary 


subsegmentation of the singular element. Twelve (12) quadrilateral subdivisions have been made 
to the singular subparametric boundary element in Fig. A. 3 and ten (10) quadrilateral 
subdivisions to the singular isoparametric boundary element in Fig. A. 4. The shaded area (0) 
corresponds to the singular portion of the element that is going to be transformed into a local 
polar coordinate system as shown in Fig. A. 5. 



Figure A.3. Subsegmentation of singular bi-linear subparametric boundary element. 



Figure A. 4. Sub-segmentation of singular bi-linear isoparametric boundary element. 



Figure A. 5. Polar coordinate transformation of singular boundary element. 


The shaded area of the singular element is fit into a new local system ( 77 , Q and subdivided into 
four triangular subelements. Each of the subelements is transformed into a new local polar 
coordinate system (p, 6) where, 

r] = pcos0 (A. 25) 

C = P sin 9 



5 * 

4 



Notice that the transformation of the differential area introduced the variable p to the integrand 
which relaxes the singularity of the fundamental solution T*(x(rj, C),0- 1° addition, an extra 
transformation is necessary to fit the limits of integration into the ( — 1,1) range such that, 


Pi = 


P2 = 


P3 = 

Pi = 


5+1 

2cos 9 
5+1 

2 sin 9 ’ 
5+1 

2 cosd ’ 
5+1 

2 sin 9 


& = £(* + 2 ) 

03 = ^(f + 4) 

04 = — (t + 6) 


(A.27) 


where the sub-indeces of the coordinates p and 6 correspond to the transformation for the 
particular integral term G 31 , G®?, G^ 3 , and G® 4 . Therefore, the influence coefficient subintegrals 
are transformed into, 

<3“ = J J r{x(i,Q,(mv,0\g^d>dt (A. 28) 

of = j'j 


which can be solved numerically by the use of Gaussian integration and through the 
corresponding transformations detailed in the previous relations. 


When dealing with isorparametric bi-linear elements, the procesure of evalaution of the element 
influence coefficients H ?• and G® ; is the same, except that the shape functions appear 

multiplying each of the integrands of the influence coefficients, see Eq. (A. 16). That is for 
instance, 

G l j = j l J 1 q t (x( V X),^)M J e (v,0\Je(v,0\dvdC (A.29) 

NG NG 

G® = y ~^W n W m q* (r) nm , Cnm, ^i)M } e (Vnm,Cnm)\Je ijlnmi Cnm)| 

n= 1 m=l 

where the subscript e is used to denote that the Jacobian is evaluated on an element basis, as in 
this case j = 1, 2...NPE = 4. Moreover, the sub-segmentation described for the G ^ is 



performed non-symmetrically depending on the location of the collocation point £* on the 
boundary element "e", see Fig. (A. 2). 
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